R/valid draft.R

In jiangrongo/HMM:

validHMM <- function(nrep, 
                     trans=matrix(c(1,1,0,0,1,1,1,0,1), 3,3, byrow = TRUE)/2,
                     M1=3){
  # M1=3  don't change. If change, also change section of violin plot
  theta <- vector("list", length=nrep)
  
  for (i in 1:nrep){
    # length of Markov chain set to multiple of M1^2+2*M1, the number of parameters
    df <- GenerateHMM(num=10, n=10*(M1^2+2*M1), trans=trans)$X
    theta[[i]] <- EstimateHMM(df, M=M1, constr=constr, tol = 0.001)
  }
  
  # initialised empty matrix for estimated means
  Est <- matrix(0, M1*nrep, 3)  
  for (i in 1:nrep){
    pos <- ((i-1)*M1+1): (i*M1)
    Est[pos,1] <- theta[[i]]$u[,  theta[[i]]$iter+1]
    Est[pos,2] <- theta[[i]]$sig[,theta[[i]]$iter+1]
    # access 3 elements in trans: 1->2, 2->3, 3->1
    Est[pos,3] <- theta[[i]]$trans[[theta[[i]]$iter+1]][c(4,8,3)]  
  }
  
  # create data frame
  df_Est <- data.frame(Est)
  colnames(df_Est) <- c("u", "sig", "trans")
  df_Est$state <- factor(rep(1:M1, nrep))
  
  # violin plot
  # to make the plot on the same level and comparable to its own true mean
  # we subtract the estimate by their true mean (Note: we set mean as state label)
  df_Est$modi_u <- df_Est$u - as.numeric(df_Est$state)     
  ggplot2::ggplot(df_Est, aes(x=state, y=modi_u, fill=state)) + 
    geom_violin() + 
    geom_abline(slope=0, intercept=0, lwd=.8, color=2)
  
  
  # plot sig2 since sig2 is less unbiased
  df_Est$sig2 <- df_Est$sig^2 
  
  ggplot2::ggplot(df_Est, aes(x=state, y=sig, fill=state)) + 
    geom_violin() +
    geom_abline(slope=0, intercept=0.1, lwd=.8, color=2)
  ggplot2::ggplot(df_Est, aes(x=state, y=sig2, fill=state)) + 
    geom_violin() + 
    geom_abline(slope=0, intercept=0.01, lwd=.8, color=2)
  
  
  # plot transition
  ggplot2::ggplot(df_Est, aes(x=state, y=trans, fill=state)) +
    geom_violin() +
    geom_abline(slope=0, intercept=0.5, lwd=.8, color=2)
  
}